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A well-posed problem in the analysis of elastic bodies 
requires the definition of appropriate constraints of the 
boundary to prevent rigid body motion. However, one is 
sometimes presented with the problem of  non-self- 
equilibrated tractions on an elastic body that will cause rigid 
body motion, while the boundary’ should remain 
unconstrained. One such case is the analysis of multi-particle 
dynamics where the solution is obtained in a quasi-static 
approach. In such cases, the motion of the particles is 
governed by the dynamic equilibrium while the contact 
forces between particles may be computed from elastostatic 
solutions. This paper presents two regularization methods of 
Interior-Constraint Boundary Element techniques for 
elastostatic analysis with improper boundary supports. In the 
proposed method rigid body modes are eliminated by 
imposing constraints on the interior of an elastic body. This 
is accomplished through simultaneously solving _ the 
governing Boundary Integral Equation and Somigliana’s 
Identity. The proposed method is examined through 
assessment and verification studies where it is demonstrated, 
that for all considered problems rigid body motion is 
successfully constrained with minimal effects on body 
deformations. 
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1. Introduction and background 


Discontinuum mechanics govern a class of problems where the materials involved are, in 
general, granular in nature or particulates, e.g., aggregates, noncohesive soils, aerosols, and 
atmospheric pollutants, among others. The numerical modeling of such problems is largely based 
on the modeling of the contact mechanics among material granules and is of common interest in 
many sciences and engineering fields. The complexity of the particle interaction in such models 
cannot be realistically represented by continuum-based models, such as the Finite Element 
Method (FEM) and the Boundary Element Method (BEM). The Discrete Element Method 
(DEM) is a discontinuum mechanics based method that finds its roots in the 1979 pioneering 
work of Cundall and Strack [1] who reported on a 2-D method for modeling granular material. 
The individual particles are customarily represented by simple shapes that move in space as rigid 
bodies. Despite the rigid particle assumption, the original form of the DEM estimates the 
particle-to-particle contact forces based on deformable particle interaction laws. The majority of 
analytical solutions that estimate the contact forces are developed from Hertzian theory. Such 
solutions exist only for a few simple shapes of the contact area (e.g. circles and triangles), and a 
review of particle interaction laws is reported in [2]. Alternatively, Mathews [3] considered 
deformable particles in a system that move and interact with each other and the system 
boundaries (if any) under the action of contact and field (non-contact) forces. The motion of the 
particles is governed by Newton’s equations of motion that are solved in a time-marching fashion 
within the DEM framework while the particle kinematic interaction forces using a quasi-static 
approach using an elastostatic BEM [3]. For the purposes of this work in problems considered 
with the proposed DEM-BEM method, bodies are set in an empty space which doesn’t contain 
any external supports. Under these conditions, if displacement or force acts on a body, the latter 
will develop rigid body modes. In the proposed DEM-BEM method an elastostatic BEM is used, 
where, it is required that all bodies are self-equilibrated and properly constrained against rigid 
body motion. In classical BEM formulations, however, the constraints are applied on the 
boundary of the solution domain. If such constraints do not stem from the physical realization of 
the problem, they tend to alter the deformation response of the body. Furthermore, it is often 
necessary that the boundary of a body remains unconstrained for the model to accommodate 
traction only boundary conditions [3]. 


The role of rigid body modes in the BEM is studied by Vable [4], where, the ways in which rigid 
body motion appears in BEM solutions is discussed and an algorithm to eliminate rigid body 
motion by modifying input data is presented. The necessary modifications are determined by 
subjecting a body to a unit rigid body motion at each boundary support and computing the 
corresponding unknown traction fields. These tractions are then used to calculate the 
displacement that needs to be “subtracted” from the input data to account for the rigid body 
motion that would have occurred in the originally considered problem. Some of the variables 
needed for this computation have to be determined from analysis for matrix conditioning 
improvement. The use of Vable’s algorithm requires several additional computations before a 
BEM system is ready to be solved. 
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Blazquez et al. [5] and Vodicka et al. [6,7] present modifications to the Boundary Integral 
Equation (BIE) to account for rigid body motion through two main approaches: a) adding 
boundary point supports (“Method S”) and b) adding integral operator(s) based on Fredholm 
Theory [8] to the BIE (“Method F’”). It is concluded that “Method S” gives acceptable results for 
problems with equilibrated loading conditions, but is not well suited for complex loading 
scenarios. Furthermore, the use of “Method S” does not allow for unconstrained boundaries. 
“Method F” approaches have been studied by several researchers as discussed in [9], and are 
based on incorporating Fredholm Integral Operators into the BIE such that invertible system 
matrices are obtained. Employment of Method F requires the use of additional variables that 
have to be selected based on several matrix invertibility conditions. 


Asadollahi and Tonon [9] use an elastostatic BEM to analyze rock block systems which have no 
naturally occurring boundary supports. Rock block systems are composed of a rock block 
surrounded on all sides by a mass of rock. The BEM is employed to compute the deformation of 
the rock block due to an all-encompassing traction field caused by the surrounding rock mass. 
Each face of the block is discretized by several triangular elements and rigid body translation is 
eliminated by fixing one arbitrary boundary node. Since only one boundary node is fixed, the 
system matrices should still be singular due to rigid body rotations. It is stated, however, that due 
to round off errors the system matrices are ill-conditioned and a matrix inverse can be found 
using the algorithm presented in [10]. This method of eliminating rigid body modes is suitable 
for the rock block systems analyzed in, [9] since, all tractions are known before the BEM solver 
is employed. This method is not applicable, however, to bodies with unconstrained boundaries. 


The removal of rigid body motion has also been examined using several other mathematical 
approaches including “Regularization” [11,12], “Singular Value Decomposition” [13], and 
“Algebraic Multigrid Methods” [14] among others. An overview of these methods is given in 
[9,15]. In general, as stated in [9], the fundamental concepts of these mathematical approach 
based methods are difficult to interpret from an engineering point of view. For problems where 
the boundary of a body should remain unconstrained, or partially constrained, while the body is 
subjected to a set of unbalanced forces, a BEM method is required where rigid body modes can 
be eliminated through internal constraints. In standard elastostatic BEM procedures, a discretized 
BIE is solved to obtain all unknown displacements and tractions on a boundary, after which, 
Somigliana’s Identity [16] is employed in a post-processing fashion to compute interior 
displacements. It is evident that interior points cannot be fixed using this solution methodology, 
and there is a relatively limited amount of literature on applying internal constraints in the BEM. 
Mina and Rashed [17] use the BEM to model floor slabs that are supported in the interior by 
columns. In this work point collocation equations are expressed for each column and solved 
simultaneously with the governing BIE. The governing equation and corresponding fundamental 
solutions used in [17], however, are based on plate bending theory which is different from 
elastostatics. In addition, the systems considered by Mina and Rashed are only loaded in one 
direction that is out of a plane with respect to the floor slabs. 


This work proposes an Interior-Constraint BEM where rigid body modes are eliminated by 
constraining points within the domain of a body. The methodology of the proposed work is 
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presented with respect to an elastostatic BEM formulation. The application to potential problems 
is a straightforward reduction on the number of constraints. It is expected that an unknown 
reactive force field, equivalent to a body force, will develop when interior points or regions are 
constrained. Subsequently, three approaches for enforcing interior constraints are presented and 
discussed in detail. The first approach assumes that the equivalent force is a concentrated force 
that develops at the constrained degree of freedom. The next approach assumes that a finite 
region is constrained and body forces develop over the finite region. The third approach assumes 
that the displacements at selected interior points are constrained, and the equivalent smooth body 
force distributions are assumed over the entire body. The proposed method is evaluated through 
assessment and verification studies where it is compared to analytical and Finite Element 
Method (FEM) solutions. 


2. Methodology 


2.1. 2D Elastostatic BEM formulation 


Consider the 2D linear elastic, homogenous, isotropic body Q bounded by Ff with outward normal 
n shown in Fig. 1. Assuming plane strain conditions and small deformations the governing 
equation is the Lamé-Navier equation expressed in indicial notation as [18]. 

: ~b; = 0 1 
(=) Uy, je + Ui jj + Oi = (1) 
where wu is the displacement field at any point, b is body force per unit mass, i, 7 = 1, 2 denote 


directions along the global Cartesian coordinate system with unit vectors e), e2, v = Poisson’s 
ratio, and “4 = Lamé’s Constant. 


2 n 2 


Tp 
e P=f,+l, 


e1 wi 


Fig. 1. 2-D Elastic Body: Nomenclature. 
The associated boundary conditions are expressed as [19], 


uj =u; onl, (3) 
tj = O{jNj = t; on lr, (4) 


where I, UT, =I’, %; and t; are prescribed displacements and tractions on boundary segments I’; 
and I, respectively, oj; is the stress field, and n; is the component in direction j of the outward 
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normal to the boundary. The associated BIE is developed through well established procedures 
[19], and is expressed as: 


chul tf. tyjujdl =f. ujtal + f, ujjbjdo (5) 


where: (i) uj and t; represent the 7 component of boundary displacement and traction fields, 


ce a” 


page (ii) b; represents the body force components, (iii) “gq” is a point on the boundary, 


66 a” 


(iv) ci is the oe Principal Value [20] at point “gq” commonly referred to as the jump term 
[19], and (v) uj;, tj 
tractions, eave The fundamental solutions are defined between any two points in the 
solution domain expressed in indicial notation as [21]: 


are the fundamental solutions (Green’s Functions) for displacements and 


ui; = SEG - In (2) 54 + Har| (6) 
tis = Gye [SE {C — 2006, + 2renj} + (A - 2v) (uy — mer,)| (7) 


where r is the distance vector between the two points, and 6;; is the Kronecker delta function 
[22]. The BIE can be expressed in a discrete form through a collocation approach [19]. The 
discretization yields N boundary degrees of freedom (DOF) and M domain DOF. The BIE Eq. 
(5) is cast then in a matrix form as: 


Hu = Gt+ Bb (8) 


Where: (i) H and G are N x N influence matrices associated with displacements, u, and tractions, 
t, on the boundary nodes, respectively, and are computed based on a collocation approach from 
integrals of the fundamental solutions over the boundary of the domain; (ii) B is an N x M 
influence matrix associated with the body forces, and is computed based on a collocation 


approach from the last integral, So uj; b;dQ, of the BIE Eq. (5); and (iti) b are the discrete nodal 


body forces. Once the BIE is ed displacements, u%‘, at points g; within the domain can be 
found using Somigliana’s identity [18] which is expressed in a matrix form as [19]: 


uti + HU = GT + Bb (9) 


where H and G are matrices of dimension M x N that contain the integration of the fundamental 
solutions over the boundary with respect to point g;, and B is a matrix of dimension M x M that 
contains the integration of the displacement fundamental solution Eq. (6), over each domain 
element with respect to point g; and involves only interior points. Details on the evaluation of the 
submatrices in Eqs (8) and (9) and the implementation of the BEM method are widely available 
in the literature, e.g. [19,22]. 


6 G.F: Mathews et al./ Journal of Soft Computing in Civil Engineering 2-2 (2018) 01-18 


2.2. Interior-constraint BEM formulation 


The proposed Interior-Constraint BEM (ICBEM) of this work eliminates rigid body motion in a 
deformable body by constraining displacements at interior points. This is achieved by 
simultaneously casting the governing BIE (Eq. (8)) and Somigliana’s Identity (Eq. (9)) in a 
matrix form as: 


a cal ft=lé Sued (10) 


Subsequently, and in view of the physics of the real problem, a number of interior points, q;, are 
selected at which the necessary displacement constraints, e.g. u7‘ = 0, should be attained. The 
problem then reduces to computing the associated equivalent body force field, b°?, that, along 
with the boundary excitations and other existing body forces, will satisfy these constraint 
conditions. Without loss of generality, it is assumed that other existing body forces are zero and, 
therefore, in Eq. (10) b = b® represents a discrete form of the equivalent body forces. It is 
evident that the discrete form depends on the spatial distribution of the body force. In this work 
three cases have been considered and investigated: (a) The first case assumes that at every 
constrained interior DOF an equivalent concentrated body force develops; (b) The second case 
assumes that a finite, yet small, area in the domain is constrained and all equivalent body forces 
are distributed over this area; and (c) an appropriate number of DOF at interior points are 
constrained and the equivalent body forces are assumed to be distributed over the entire domain. 
As discussed in the following sections the case of concentrated body force leads to strong 
singularities and a solution cannot be obtained, while the other two cases lead to stable and 
accurate solutions where equilibrium is satisfied. 


2.2.1. Concentrated equivalent body forces 

The case of concentrated body forces (i.e., body force distribution is a delta function) is 
considered first to demonstrate that this option leads to singular integral values. When a single 
interior point is constrained the body force term of the BIE (Eq. (5)) becomes unbounded. This is 
demonstrated by considering the computation of the body force term with respect to the circular 
area shown in Fig. | when loaded by a constant body force, i.e. 


J, UijbjdO = 


ae ah i |G- 4v) In (=) 5; + rary] dédr (11) 


Where r and @ are the radial and angular coordinates and the resultant body force (R,) acting 
over the circle is expressed in terms of the constant F' as 


E 2 F 
R, =f, fy 3 40dr =F (12) 


When ¢ > 0, Eq. (11) is equivalent to computing the body force term with respect to a single 
point, i.e. 

‘ 2 1 
J, ujibjdO = limeo {Jy fo [C3 — 4v) In (=) 54 + rary |= adr} (13) 


1672 na v) 
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The derivative terms in Eq. (13) are non-singular and can be evaluated without any issues. 
Evaluating the other terms in Eq. (13) the following expression is obtained when i =: 


wee lime-so { a [In (-) 5; i| ~ ddr} = es limeo {1 +In (=)} (14) 


It is evident that Eq. (14) does not converge to a finite value when ¢€ approaches zero. The same 
analysis of computing the body force term over a single point has also been considered for linear, 
polynomial, and exponential body force distributions. In all cases obtaining a non-singular 
solution is not feasible. Special Gauss Integration schemes (e.g. [23]) which are often employed 
in the BEM to treat singular integration over boundary elements cannot be used to evaluate the 
integrals in Eq. (14), since, the area of integration tends to zero. There are some mathematically 
based approaches for approximating divergent integrals by eliminating terms that tend towards 
infinity such as Hadamard Regularization [24], but this is beyond the scope of this work since it 
increases the complexity of the computations. It is therefore concluded that in the ICBEM all 
body forces need to be distributed over a finite area, as discussed in the following two sections. 


2.2.2. Body forces distributed over a small, finite area 

This approach assumes that a finite, yet small, area in the domain is constrained and all 
equivalent body forces are distributed over this area. In this work, the constrained area is a 
Triangular Region (TR), as shown in Fig. 2. The displacements at the TR vertices are 
constrained, u4i = 0, and the boundary is discretized with quadratic line elements. The discrete 
body force vector b in Eq. (10) consists of values of the body force at the triangle vertices. The 
TR is sized in such a way that its effect on body stiffness is negligible, and this is demonstrated 
in the Assessment and Verification Studies section. 


Boundary Discretization 


@ End Nodes O Middle Nodes 


Fig. 2. ICBEM discretization for body forces distributed over a finite region. 


2.2.3. Body forces distributed over the entire domain 

In this approach, in order to remove the two translational and one rotational rigid body modes, 
two interior points are defined first. The first interior point, point “c”’, is typically placed, without 
loss of generality, at the body’s center of gravity and the second point, point “o”, is placed at an 
offset distance from point “c” in any direction, as shown in Fig 3a. The displacements u° = 


66 99 


{uf uS}" at point “c’”, shown in Fig. 3b, are set to zero to constrain rigid body translation. The 
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associated body force field, b‘’ = {b{" b5"}", is assumed uniformly distributed over the entire 
domain as shown in Fig. 3a, where, “tr” denotes that the body forces are associated with the 
translational constraints at point “c”. Rigid body rotation is constrained by restraining the 
displacement at an additional interior point “‘o” in a direction perpendicular to a line that passes 
through “c” and “o, as shown in Fig 3c”. In the implementation of the method in this work, and 
without loss of generality, point “o” is offset horizontally from point “c”. Consequently, to 
constrain rigid body rotation the vertical displacement, wu? at point “o” must be constrained. To 
develop the ensuing internal moment from constraining rigid body rotation, a rotational body 
force field, b’,with a linear distribution over the entire domain about point “c” is assumed, as 


shown in Fig. 3a. 


(b) (c) 


Fig. 3. ICBEM methodology for body forces distributed over entire domain: (a) Body force fields, (b) 
Translational constraints, (c) Rotational constraint. 


Introducing the displacements wuj,us, andu? and the associated body force fields 
bi", bs’, and b" in the ICBEM Eq. (10), the latter becomes: 


H —B —-B’](U G 0 0 U 
Ho —Btre —Brcl pts = G° | 0 ut (15) 
H° —Btro —Bpro bt G° 0 =A us 
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where: (i) the first equations are the discrete BIE, (ii) the middle equations are Somigliana’s 
Identity for the constrained displacements, u°, at point “c”, and (iii) the last equation is 
Somigliana’s Identity for the constrained displacement, u9, at point “o”. Matrices G°, H°, G° and 
H°, are the H and G matrices defined in Eq. (10) and evaluated with respect to internal points “c” 
aaa “9”, respectively. Matrices B™*° and B'™° are of size 2 x 2 and 2 x 1, respectively, and 
matrices B™* and B® are of size 1 x 2 and 2 x 1, respectively and are defined based on the body 
force term in Eq. (5) where the displacement fundamental solution is evaluated with respect to 


66 99 


points “c” and “o”, as discussed next. 


2.2.3.1. Numerical evaluation 

The evaluation of these terms requires integration over the domain of interest; however, a 
discretization of the domain is not necessary. Instead, the domain is divided into triangular 
regions, as shown in Fig. 3b and 3c that represent convenient, simple areas over which Gaussian 
quadrature can be applied. Each triangular region is defined by point “c” and two consecutive 
boundary nodes. In general, the most suitable physical representation of a problem is achieved by 
placing point “c” at a body’s center of gravity, and this is done for all problems in this work. In 
addition, for this work, each triangular region is defined by one integration point located at the 
centroid of the region area. 


In view of the assumption that the equivalent body forces associated to translational constraints 
are uniformly distributed, the elements Bi in matrices B™° and B'™° are evaluated numerically 
as: 


Jig vidya = EET LE ui dO} DE = BEDS (16) 


where Q, is the area of a triangular region. For the one point, Gauss quadrature employed in this 


work, Eq. (16) the Bj; Bit terms are evaluated as: 


Br= = yeas THOR (17) 


Where uj; is shown in Eq. (5) for a distance vector defined between points: (i) “c’” and the 


are) bed 


centroid of each triangle for the B™* termthe, and (ii) points “o 
for the B™° term. 


and the centroid of each triangle 


In view of the linear distribution of the equivalent rotational body force, b’, the latter is written 
in vector form as 


b’ = —[b’esinale, + [b’e@ cos ale, (18) 


Where: (1) b” is the magnitude of the rotational body force field, and (ii) a is the angle (measured 
counter-clockwise) between the horizontal, and a distance vector with magnitude @ that is 
defined from point “c” to the integration point of each triangular region (depicted in Fig. 3c). Ina 
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manner similar to the evaluation of matrices B™’* and B™°, the elements B/ in matrices B™, and 
B”° can now be expressed in terms of one point Gauss quadrature as 


pr —_ wt Triangles 
B; ~~ pee 


(—uj, sina + uj, cos a)eN, (19) 
Finally, matrix B” accounts for the effects of the rotational body force field in the BIE and is 
computed in the same way as Eqs. (19), but the source points correspond to boundary nodes 
instead of the constrained domain points “c” and “‘o”. The BEM system matrix is square with 
dimension N+3 x N+3. The solution to Eq. (15) will yield all unknown field variables, after 
which, the equilibrium of a body can be verified with the following equations 


bf’ QO+F, =0 (20) 
ya Triangles pr 20+ M =0 (21) 


Where F’, and F> is the resultant external force on the body in the 1, 2 directions, respectively, 
and & is the total external moment applied on the body about an axis that is orthogonal to the 1, 
2 directions. 


3. Assessment and verification studies 


To assess the proposed ICBEM method, it is examined through showcase problems and 
parametric studies where it is compared to other solution methods. In this section the two 
variations of the ICBEM will be denoted as methods: (1) M1 — Rigid body motion is constrained 
by fixing a single domain region, TR as discussed in section 2.2.2, and (ii) M2 — Rigid body 
motion is constrained by fixing a combination of degrees of freedom at two domain points, and 
body forces are distributed over the entire domain (section 2.2.3). The following studies are 
performed: 


(i) Methods M1, M2 are employed to solve the problem of a steel plate symmetrically 
deformed in tension. In this study, the ICBEM is compared to a known analytical 
solution and the FEM. 

(ii) Methods M1, M2 are further studied by reconsidering the steel plate problem with an 
asymmetrically applied displacement field, where, no analytical solution exists for this 
set of boundary conditions. 

(iii) A parametric study is conducted to assess the effects of the size of the constrained 
region, TR, with method M1. This study is performed considering the steel plate 
problem in the preceding study. 

(iv) Methods M1, M2 are employed to solve the steel plate problem when loaded by a 
traction field which produces an external moment. The accuracy of the ICBEM is 
assessed by examining plate deformations and through force and moment equilibrium. 
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3.1. Steel plate deformed in tension 


The 2m x Im steel plate shown in Fig. 4 is deformed in tension on both sides by u = 0.001m. The 
plate material has a modulus of elasticity E = 200,000 MPa, and Poisson’s ratio v = 0.25. The 
boundary mesh used for methods M1, M2 is shown in Fig. 4(a), where the perimeter is 
discretized in 120 quadratic boundary elements. The triangle at the center of the plate represents 
the fixed TR area employed in method M1. For this initial study, the area of the triangle is A, = 
0.0015 m* which corresponds to 1/1333.33 of the area of the plate. In method M2 for domain 
integration, there are 2 triangular regions for each quadratic boundary element which 
corresponds to 240 triangular regions in total. Point “c” in method M2 is placed at the center of 
the plate, and point “‘o” is offset 0.05 m to the right of point “c”. The problem is also solved by a 
plain strain FEM analysis (performed with ABAQUS [25]) that is designed to simulate the 
physics of the model solved with method M1. To this end, a TR area is fixed at the center of the 
plate in the FEM analysis, and the nodes of this element are denoted by the yellow circles of the 
FEM mesh shown in Fig. 4(b). The solution obtained from this FEM analysis is denoted as 
FEM1. It can be seen that even for a simple geometry a much finer mesh is required in the FEM 
than with the ICBEM. The deformation of the plate (scaled up by 1000) from the ICBEM and 
FEM solutions are shown in Fig. 4(c). It can be seen that the deformation is almost identical 
between methods M1, M2 and FEM1. In addition the deformation matches the physical behavior 
of a plate pulled in tension. It is observed that no significant deformations are introduced by the 
presence of a fixed TR area in method M1. 


> 
> 
> 
A u > u 
Fixed TR ce 
> 
> 
> 
(a) (b) 
Undeformed Plate Deformed Plate (x1000) 


(c) 
Fig. 4. Steel plate deformed in uniaxial tension: (a) ICBEM Mesh, (b) FEM Mesh (c) Deformed 
configurations as obtained by M1, M2, and FEM solutions. 


The magnitude of the resultant horizontal force on the left side of the plate has been computed 
with all methods and is reported as Case / in Table 1. The force in Table 1 is compared to the 
theoretical solution for plain strain conditions [26], 1.e. 
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Fp = 0,4p = [1 -v) - Send ae 
pi Sie GaGa) be ea atc (22) 
where 0;, €, are the stress and strain, respectively, in the horizontal direction on the left side of 
the plate, and Ap denotes the plate cross sectional area. The percent difference between all 
methods is negligible at less than 1%. 


Table 1 
End Force Comparison Study. 
Method F,, (MN) 
M1 212.91 
FEM1 213.81 
M2 212.52 
Eq. 22 213.33 


To verify that rigid body motion has been constrained by the ICBEM, displacements at the fixed 
domain points have been back-calculated using Somigliana’s Identity for methods M1, M2. 
These computed displacements are on the order of 1 x 10°!” = 0, and this indicates that the plate 
has been successfully constrained against rigid body motion. In this example, an arbitrarily small 
value was selected for A, in method M1, but almost identical results are obtained using much 
smaller TR areas. For example, using A, = 1.5x10° m? yields Fp, = 212.54 MN. It will be 
observed in the forthcoming sections that this is only the case for systems with equilibrated 
loading conditions. For method M2 the same value of F,; is obtained when point “o” is offset in 
a variable direction and distance from point “‘c”. The offset of point “o” does however affect the 
deformation behavior of the plate. To avoid introducing unwanted deformations in method M2, it 
is recommended that point “‘o” be placed as close to point c as possible. 


3.2. Steel plate with an asymmetrically applied displacement field 


The steel plate problem of the previous section is reconsidered with the applied displacement 
field u = 0.001m applied on the left side of the plate only and with no other boundary constraints 
as shown in Fig 5(a). The problem is solved with the ICBEM and FEM methods. In this example 
two FEM solutions denoted by FEM1, FEM2 are presented, and these solutions correspond to 
the emulation of methods M1, M2, respectively. For methods M1, M2 and FEM] the problem 
setup is the same as in the last section with the exception of the new boundary conditions. In 
order to obtain solution FEM2 that is consistent with the model that corresponds to M2, it is 
necessary to load and constrain the plate as shown in Fig. 5(b). In this case, the plate is loaded 
with the body force field b™ obtained from the method M2 solution, and the left end of the plate 
is constrained against horizontal displacement. In addition, the mid-point of the left end is 
constrained against vertical displacement to ensure that all rigid body modes are accounted for. 
Figure 5(c) shows the plate deformation, where it can be seen that the physical behavior of 
methods M1 and FEM] are significantly different than that of methods M2 and FEM2. 
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Undeformed Plate 


EEEEED) 


(a) 


b™ from M2 solution 


(c) 
(b) 
Fig. 5. Steel plate in asymmetric deformation: (a) Induced deformation for M1, M2, and FEM1 solutions; 
(b) Boundary conditions for FEM2 model only, (c) Deformed configurations as obtained by all methods. 


As expected the deformations from methods M1 and M2 are almost identical to their counterpart 
solutions FEM1 and FEM2, respectively. In methods M1 and FEMI the deformation is non- 
uniform, whereas, for methods M2, FEM2 the plate deforms in a uniform tapered manner. For 
methods M1, FEM1 the deformation is localized about the constrained domain region over 
which all body forces are distributed. For methods M2, FEM2 the plate deformation is uniformly 
tapered in the horizontal direction which is analogous to the distribution of the uniform body 
force field b' over the entire domain. The deformations obtained from both ICBEM 
formulations are accurate within the assumption of how body forces are distributed over the 
domain. The magnitude of the resultant horizontal force on the left side of the plate has been 
computed for all methods and is reported in Table 2. 


Table 2 
End Force Comparison Study. 
Method F,, (MN) 
M1 140.41 
FEM1 138.15 
M2 291.80 
FEM2 289.10 


It should be noted that this force cannot be computed analytically with respect to the loading 
condition of the plate. As expected the difference in force is negligible between methods M1, M2 
and their counterpart finite element solutions FEM1, FEM2, respectively. There is a significant 
difference in force between methods M1 and M2, and this is attributed to the difference in 
assumed body force distributions of the two approaches. The force from method M2 is notably 
larger than in method M1, since, in method M2 the body forces are distributed over a 
significantly larger area than in method M1. The displacements at the fixed domain points in 
methods M1, M2 have been back-calculated and again are on the order of 1 x 10° = 0 indicating 
that rigid body motion has been eliminated. Using a notably smaller triangle size of A, = 1.5x10° 
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° m? in method M1 yields Fp, = 85.43 MN. It can be seen that this force does not match the 


values reported in Table 2, and a correlation between Ay and F,, is evident. This correlation is 
studied in the next section. 


3.3. Method M1 — Parametric Study 


To further study the ICBEM the selection of A, in method M1, is examined through a parametric 
study. The steel plate problem with the unbalanced load shown in Fig 5(a) in the previous section 
3.2 is chosen for the parametric study since the resultant force on the fixed TR will be different 
than zero. The plate boundary is discretized into 120 boundary elements, however, the area of the 
fixed TR varies from A, = 1x10“ m? — 0.02 m’, representing 0.005 % - 1% of the plate area. For 
comparison, the finite element solution FEM1 (see sections 3.1 — 3.2) is also computed for all 
Ay. In all solutions, the magnitude of the horizontal component of the resultant force, F,,, on the 
left side of the plate is computed. In addition, a FEM solution is obtained for a reference 
problem. In this case, the plate is fixed at a single point at its center against translations. The 
reference problem represents the limiting case where the constrained area TR approaches zero in 
the M1 and FEM] solutions. The horizontal component of the resultant force on the left side, F,, 
is determined and assumed the reference solution. The findings are summarized in a graph form 
that relates the normalized force F,,/F, to the normalized area Ap/A, and is shown in Fig 6. It is 
observed that M1 and FEM] solutions agree well for the all areas considered. It is evident that 
the convergence to the reference solution is inversely proportional to the size of A,, since the 
rigid TR area affects the stiffness of the plate. Larger areas introduce errors while extremely 
small areas may lead to instabilities due to roundoff errors. As a rule of thumb the size of A, 
should be selected in the range 0.001 % to 0.01% of the area of the domain to ensure that: (i) The 
rigidity of the fixed TR area does not affect the stiffness of the domain and (ii) the fixed TR area 
is large enough that its effect can be numerically accounted for in method M1. 


2.0 
1.8 O FEM1 


1.6 — Mi 
1.4 


Fai Fe 


1.2 
1.0 
0.8 
0.6 


0.4 
O 025 05-075 12 125. 15 17° 2 
AD/Ay x 10000 
Fig. 6. Method M1 — Effects of fixed TR size on force solution. 
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3.4. Steel plate loaded with traction field that produces a moment 


For this study, the steel plate problem is reconsidered when loaded by the traction field t = 324 
MN/m shown in Fig. 7(a). This problem is solved with methods M1 and M2. The deformation of 
the plate for both methods is shown in Fig. 7(b)-(d). For method M1 when Ap/Ayj is selected 
within the minimum and maximum limits established in the last section the plate deforms as 
shown in Fig. 7(b). It can be seen that this deformed configuration resembles a rigid body 
rotation and is not representative of the loading condition. The plate is in force equilibrium, and 
the displacements at the fixed TR area are on the order of | x 10°'° indicating that rigid body 
motion has been constrained. It appears, however, that the plate does not properly develop the 
external moment produced by the traction field for a selection of Ap/A, within the established 
limits. To help develop the moment the problem is reconsidered using method M1 with two fixed 
TR areas. The additional fixed TR area is offset to the right of the original fixed TR area (at the 
center of the plate) by 0.2 m. The corresponding deformed configuration is shown in Fig. 7(c) 
where it is compared to a solution from method M2 shown in Fig. 7(d). The deformations are 
fairly similar in behavior but differ in magnitude. The localized deformations in the plate due to 
the traction field can clearly be seen in method M2. It can be seen that the deformation is 
representative of the applied traction field, and the plate undergoes both translational and 
rotational deformations which indicates that no rigid body motion occurred. For method M1 the 
deformation still largely resembles a rigid body rotation. It is concluded that method M1 is better 
suited for applications where rigid body rotations are expected to remain small. In method M2 
the external moment is balanced by the rotational body force field b" (Eq. (21)). 
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Fig. 7. Deformation of steel plate with applied traction field: (a) Plate loading; (b) Method M1 for 
Ap/A, between minimum and maximum limits; (c) Method M1 with two TR; (d) Method M2. 
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Force and moment equilibrium for the plate have been verified using Eqs. (20) — (21) and the 
results are reported in Table 3. It can be seen that the internal forces and moment obtained from 
the integration of body forces are equal and opposite to the external forces and moment obtained 
from the integration of tractions. The residual forces and moments are negligible, as shown in the 
last row of the table, verifying, thus, that the body is in equilibrium. 


Table 3 
Verification of Plate Equilibrium. 
F, F, M 
(MN) (MN) (MN-m) 


External Forces 
Traction Integration 
Internal Forces 
Body Force Integration 


Residual Forces 0.048 0.002 -0.040 


-167.400 329.400 206.550 


167.448 -329.398 -206.59 


4. Conclusions 


This work discussed in detail the development of an Interior-Constraint BEM (ICBEM) for 
regularization of elastostatic problems with improper boundary conditions, which in general, 
exhibit rigid body motion. The regularization process generates an unknown reactive force field 
due to the imposed interior constraints that is equivalent to a body force field. Three assumptions 
related to the equivalent force field are investigated: (i) a concentrated force develops at a 
constrained degree of freedom; (ii) the force field is uniformly distributed over a finite, yet small, 
triangular constrained interior region, method M1; (iii) the displacements at interior degrees of 
freedom are constrained and the equivalent body forces are distributed over the entire body, 
method M2. The evaluation of the proposed approaches has shown that: 


a. The assumption of concentrated body forces at discrete constrained degrees of freedom 
leads to integration singularities and is not valid. 

b. Methods M1 and M2 successfully remove rigid body translations and satisfy static 
equilibrium. 

c. Both methods have been successfully verified through comparisons with independent 

solutions based on the Finite Element Methods. 

Methods M1 and M2 correspond to solutions of different physical systems. 

Method M1 is simpler to implement than method M2. 

f. Method M1 is better suited for applications with rigid body translation only, when a 
single interior region is constrained. For applications with rigid body rotations a separate 
interior region needs to be constrained. 

g. The relative size of the constrained interior region in method M1 to the size of the 
solution domain affects the solution. These effects are more pronounced in the 
computation of domain forces and less evident in the deformations. Recommendations 
on the relative size have been established. 

h. Method M2 is shown to be more robust and independent of solution parameters. 


a 
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